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Abstract 



A kink-based expression for the canonical partition function is developed 
using Feynman's path integral formulation of quantum mechanics and a dis- 
crete basis set. The approach is exact for a complete set of states. The 
method is tested on the 3x3 Hubbard model and overcomes the sign problem 
seen in traditional path integral studies of fermion systems. Kinks correspond 
to transitions between different N-electron states, much in the same manner 
as occurs in configuration interaction calculations in standard ab initio meth- 
ods. The different N-electron states are updated, based on which states occur 
frequently during a Monte Carlo simulation, giving better estimates of the 
true eigenstates of the Hamiltonian. 
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I. Introduction 



Studies of disordered and/or finite size electronic systems, sucli as clus- 
ters, amorphous solids, or quantum dots pose challenges for computational 
methods. The difficulties include the need to accurately treat electron- 
electron correlation and to include finite temperature effects (particularly 
as they affect the atomic positions.) The present work is motivated by the 
desire to simulate large, multi-electron clusters. Studies of these systems are 
hindered by the need to identify and/or explore global and local minima. 
The rapid growth in the number of minima with cluster size makes the de- 
velopment of precise, accurate and fast computational algorithms essential 
for the study of large systems. Often, isomers have similar energies, which 
requires both the calculation of free energies and the accurate inclusion of 
electron-electron correlation in order to have a viable method. Feynman's 
path integral (PI) formulation of quantum mechanics has the potential to 
provide such a method and is the focus of this work, as its features include 
exact inclusion of correlation and the calculation of the partition function, 
which allows both the correct sampling of different geometries and the si- 
multaneous treatment of electronic and geometric degrees of freedom. The 
latter is a major advantage, as it does not require an accurate estimate of the 
electronic energy at each set of atomic coordinates. To see this, we can write 
the partition function for a system with electronic and geometric degrees of 
freedom as 



where R'^ denotes the set of geometric coordinates, {a} denotes the elec- 
tronic basis set, and p is the canonical partition function at an inverse tem- 
perature f3. If this is evaluated by sampling {R^} and {a} from \p\, we see 
that new geometries can be sampled without having a converged electronic 
state {{a} does not need to correspond to the ground or excited state). 
Effectively, the electronic and geometric optimizations can be carried out 
simultaneously, in contrast to conventional electronic optimization methods, 
in which a converged energy or set of forces is necessary when changing ge- 
ometries. The result is that, in principle, the sampling of geometric phase 
space can be accomplished in a time on the order of the time required for a 
single electronic energy calculation. 



Q 




{a} 



(1) 
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Most path integral approaches use a complete set of states, namely, the 
position representation, which offers a route to an "exact" determination of 
the partition function for a many-body system.Biiiii'ii'illl3ii'0ili0i 
These approaches have been frustrated by the sign of the fermion density 
matrix, which can be positive or negative and can lead to computational 
difficulties in obtaining precise results. Thus, before the PI method can be 
applied to the study of clusters and other multiple-minima systems, it is 
necessary to develop ways to avoid the sign problem. The sign problem is 
independent of whether or not the atomic positions are fixed. Hence, in this 
first paper on the subject, we focus on electronic degrees of freedom with 
fixed atomic positions and leave the extension to varying atomic positions 
to future work. We note in passing, though, that many interesting systems 
have fixed atomic positions (e.g., quantum dotsEl), so that the developments 
in this paper can have immediate applications. 

Herein, we abandon the use of the position representation and use a finite 
basis set in order to avoid some of the difficulty inherent with the position 
representation. It is well known that finite basis sets (for example, Gaussian 
basis sets) are capable of producing accurate results for many systems. § We 
use finite basis sets and the discretized version of the path integral expression 
for Q, the canonical partition function, to develop a "kink" expansionillii for 
Q. In the discretized version of Q, paths are divided into small imaginary time 
segments. When using finite basis sets a path will spend some imaginary time 
in one (many-electron) state (this time may involve several imaginary time 
segments), have a transition to another state during one imaginary time seg- 
ment, spend some imaginary time in the second state, have a transition, etc. 
The transitions between states are called "kinks" .il A path can therefore 
be classified by the number of kinks and states involved. When we rewrite 
the expression for Q in terms of kinks and states, we call this a kink ex- 
pansion. We note that kinks can correspond to excitations from the ground 
(Hartree-Fock) state analogous to those seen in configuration interaction (CI) 
calculations. The expression we develop will give the exact value of Q (in- 
cluding electron-electron correlation) for a complete set of states. The zero 
kink contribution to Q will have no sign problems, since the system is in a 
single state. With a properly chosen ground state, it is possible that conver- 
gence of the path integral can be obtained with just a few kinks; this would 
significantly reduce the sign problem, thereby increasing the speed and pre- 
cision of a calculation. A good choice for the ground state can be obtained 
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using a Monte Carlo simulation, in which the different N-electron states that 
appear during the simulation are used to update the estimates of the ground 
and excited states in a way that includes electron-electron correlation. We 
call this approach an adaptive approach, since the Monte Carlo algorithm 
allows the estimates for the ground and excited state wavefunctions to evolve 
according to the statistical sampling of the different N-electron states. 

In the following, we develop a closed form expression for the canonical 
partition function, cast as a kink expansion. We apply our formulation to a 
model problem, the 2-D Hubbard model and demonstrate the efficacy of the 
approach. 



II. Kink Formalism 



The kink formalism described here assumes a set (usually finite) of or- 
thonormal, N-particle states, which we denote by {oi}. In terms of these 
states, the partition function can be written 

Q = Tr{e-/^^} = E(«.|e-^^|«,) (2) 

3 

We write this as 

Q = \mvQ (P) 



P^oo 



Q{P) = ("ii|exp(--i/)|ai2 ) ("i2|exp(--//)|aj, 




/9 



(ajp\e^v{-^H)\aj}^ (3) 

The introduction of P allows high temperature approximations for the matrix 
elements. In this form, the sign problem can be easily seen. Since each of 
the matrix elements in Eqn ^ can be positive, negative, or zero, the sign 
of the summand can be negative for some sets of {aj}- Thus, during a 
Monte Carlo simulation, the estimator for Q can change sign, leading to 
large statistical errors. This alternation of sign is the source of difficulties 
in using the path integral formulation to study multi-fermion systems. We 
focus on an expression for Q (P). The kink expansion is obtained by recasting 
the sum over {aj} as a sum with all jk equal (no kinks), one jk different (2 
kinks), two jk different (3 kinks), etc: 
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Q{P) = E 



aj^\exp{-^H)\aj^ 



+ 



P-2 

EE 

ii J2 n=0 



aj^\ex.p{-^H)\aj, 



«j2|exp(--// 



\a 



32 



P-2-n 



P 



1 2 



«ji|exp(--//)|aj2 



+ 



X 



(4) 



where the first term is the zero kink term and the second the two kink term. 
In this and the following equations, ji ^ j2, J2 7^ js, etc. If we set 



e = P/P 



(5) 



Xj = {aj \ exp{—eH)\aj 



and 



(6) 
(7) 



tij = {ai\ exp{-eH)\aj) 
we can write the expression for Q (P): 

Pin \ / " \ I P^n\ / « \ 

l+l2 + -+ln,P- 

j n=2 \i=l ji ) \k=\ / \k=llk=0/ \k=l / 

(8) 

where the first term is the zero kink term and jp+i = ji- Recognizing that 
there are (^^jways to put the n kinks at the different P sites, we can choose 
the location of the first kink and rewrite our expression as 

Q{P) = E^f + 



E^ HE n^... 



n=2 \i=l ji I \k=\ 
P-nP-n-l„ P-n-ln-ln-i ^3 

E E ••• E 



X 



3n Jn-l J2 31 



'(9) 



i„=0 i„-i=0 



^2=0 
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With the shorthand notation Sj = ln + ln-i-\ \- Ij: we have 



P p / n \ / 

E- HE 

n=2 \i=l ji ) \k=l 
P-n P-n-S„ P-n-S:i 

E E ■■■ E ^"-"^ 

ln=0 l„-l=0 l2=0 



X 



lnJri-1 _ _ _ h P-n- 



Consider the term 



P I " 



p-n P-n-5„ P-n-Sz 



n 



=1 3i I \k=l J Z„=0 «„_i=0 /2=0 



Qn = - HE nw E E ••• E 



We first assume that Xj^ ^32 ^ ' ' ' ^jn define 

P-n P-n-Si 



Using, 



we find 



S{{xj},n) = J2 ■ 

lrt=0 



P-n-S3 



"32-"3l 



l2=0 



J2 xi-x""' 



J2 ^P-n-S2 

3n'"3n-l •^32-^h 



l2=0 



M+1 



P-n-S.i 

E 



P-n-53 



i2=o 



31 



X 



p-n-53 



P-n~S:i+l- 



P-n-53+1 _ „P-n-53+l 

2 P-ri-S:t + l 

E^-^^- — 

n {xji, Xj^) 

S{{x,},2,n) 



We now proceed by induction to develop a general form for S {{xj} ,n). 
Assume that 

i_l P-n-5.+{i-2) 

S{{x,},^-l,n) = J:-Jf 



Consider the next summation in Eq. |l^: 

P-n-Si+i i-1 ^P-n-Si + {i-2) 



Now 



n {x,,-x,^: 



P-n-Si+i 



H XuS{{xj} ,i-l,n) 

P-n-Si+i i-1 ^P-n-Si+i-li + (i-2) 

E 



=0 U{x,,.-x 



i-1 

I nr 



i-1 P-n,-5,+i + (»-2) 1 _ 
2^ i-1 ^ ^ _ 

kj^m 

i-1 P-n-Si+i + (i-l) i-1 i-2 P-n-Si+i+1 

E^f E — ^"^-^ (15) 



^ n {xjf, Xj^) ^ 1 (xjj, Xj-) n {x 



k^m k^m 



i-1 i-2 P-n-Si+i+1 



i-1 

t nr . — T . 



i-1 

^—^ {Xji_ — Xj-) n {Xji_ — Xj^) 

P-n-5,+i + (i-l) i-1 ^ UiXj^-XjJ 

n (^ji ~ Xjf,) ^-^ (xj^ — Xj^) n {xj^ — Xj^^ 

kjt^i niytk 
P-n-5,+i + (i-l) i-1 Qf-Y ^ n (Xj^-XjJ 

xE- 



i-1 



u{xj,-x,,) ^==1 n(x, 

fc^^i m^k 



3k 
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i-1 



P-n-Si+i+{i-l) i_i n (Xj, - Xj^ 



-^E=^ (16) 



Notice that 



i-i n {xji — Xj^) 

E^^ -E^.(^.J (17) 

^==1 n (1 - ^) 

i-l 

is an i-2 order polynomial in Xj-. Consider y{x) = Yl h{x). This is an i-2 

fc=i 

order polynomial in x. Since lk{xj^) = x'^~'^6k,i, y{x) is an i-2 order polynomial 
that has the value at the i-l points Xj-^^,Xj^, . . . ■,Xj^_^. Since y{x) is an 

i-l 

i-2 order polynomial, we must have lk{x) — x . Therefore, 

k=l 

i-l n {xji — Xj^) 

Et^ = 1 (18) 



i-2 i-l 

Xji k=l Yl {1 — — 
m^k 



and 



i-l P-ri-'^^ni + a-^) i-l ,i-2 P-n-S,n + l 

" ^ ' Jk ' J I 



E i E i-l 

'^-i n {xj^ — Xj^) [xj^ — Xjj n {xj^ — Xj^) 

ky^m k^m 
i P-n-Si+i+{i-l) 

= E^ (19) 

^-^ n (^^ifc ~ Xj^) 

k^m 

= 5 {{xj} , i, n) 
Therefore, by induction, we can show that 

S{{xj},n)^±^^^ (20) 

fe=i n {xj^ - Xjj 

k^m 



9 



Next consider the case where some of the Xj^, are equal. The general result 
can be inferred by assuming Xq occurs twice in the sum: 

P-m P-m-k P-m-k-Si P~m~k-Sn-i P~m-k-Sn P-m-k-Sn-j 

■ ■ ■ E E E ^2 ■ ■ ■ E E ^0 E 

k=0 h=Q l2=0 /„=0 j=0 1=0 

(21) 

M M-l M M-k 

E^iE4 = E^oE^I (22) 



Now 



/=o fe=0 fc=0 /=o 



so Eq. becomes 



P-m P—m—k P-m—k—j 

■■■E^o E 4 E ^i'--- (23) 

k=0 j=0 h=0 

As a result, if there are s identical xo, we will arrive at 

p-m P-m-Si P-m-S2 P-m-Sk^-i P-m-Sk^ 

••■EE E ••• E xo^^+^^+-+^^ E ^r--m 

ki=0 k2=0 fc3=0 ks=0 1=0 

P-m P-m- Si P-m- S2 P-m-Sks-2 P-m P-m-ks 

= •■•E E E •■■ E E 4^ E -!■•■ 

ki=0 k2=0 fe3=0 ks-i=0 ks=Sks-i 1=0 

P-m ks ks—Si t^s—Skg—2 P—m—ks 

= •••E^^E E ••• E E ^\--- 

ks=0 ki=0 k2=0 ks~i=0 1=0 

where we used Ek=o^ EfLk+N = T.f=NT.i'JS ■ Let 

ks ks—Sl ks—Sks-2 

w{s,h)^Y: E •■■ E (25) 

ki=0 k2=0 ks-i=0 

Now we assert that W{s,k,) = [''''^'^^). When s = 2, W{s,k,) = fc, + 
1 = (^'^s'^]^^)- To show this in general, we use proof by induction. If 
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^ (fc, + g-2-fci) 
fc'^o(^-2)!(A;.-A:i): 
(g-2-A;i)! 
(--2)! MO! 

^ (g-2 + A;i)! 
,^J.-2)!(A:0! 



(26) 



where the last identity is taken from Gradshteyn and Ryzhik.il Thus, Eqn. 
E3| becomes 



P-m / h. \ a l\ P-m-ks 



ks=0 



1=0 



P-m ^s-l 



P—m—ks 



„ks+s—l 



E 

«=0 



a;, 



(27) 



Thus, we have for m distinct x^^'s, each occurring Sj^, times 

m 

S ({xj},n,m,{sj}) = n 



k=l 



1 ci^Jfc-^ 



(s,, -Olrfx^- 



m ^P—n+m—l 
1=1 n {^ji ~ ^ji) 



The derivatives can be evaluated recursively and therefore we have a final 
expression for Q: 



QiP) = Y.^f + 

j 

fnEl fn S{{x,},n,m,{s,}) (29) 

n=2 \i=l ji J \k=l / 

This expression can be evaluated using a Monte Carlo algorithm in which 
both N-electron states and kinks are sampled. While the number of kinks 
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can grow to be P, in practical calculations it is plausible that a good choice 
of states will require only a limited number of kinks. If there are or 2 kinks, 
it is clear that the Monte Carlo estimator will be positive; for more than 2 
kinks, there will be some negative values for the estimator. However, as the 
estimate for the ground and excited states improves, the sign problems will 
be reduced. 

III. Application: 2-D Hubbard Model 

The formalism was applied to the 2-D Hubbard model.@ In this model, 
the Hamiltonian is given by 



where t is a hopping term that allows hopping between nearest neighbor 
sites and U is an repulsive energy term that has non-zero contributions when 
two particles are on the same site. In this work, we chose a 3x3 lattice, 
a = — 1.0, 6 = 4.0, and periodic boundary conditions were not used. We 
applied the kink formalism to several occupancies of the lattice: Neiectron = 
3, 4, and 5. The number of up- and down-spin electrons is given in Table |l|. 
Temperature was chosen to represent "typical" temperatures encountered in 
studies of clusters: jSAE > 10 where AE is the difference in energy between 
the ground and first excited state. This value was chosen assuming a 0.1-1 
eV gap between ground and excited states (a typical value in small metallic 
clusters) and a temperatures on the order of a hundred Kelvin). The values 
of /? are given in Table |1|. 

The method was made adaptive in the following manner. An initial set of 
N-electron states was chosen by diagonalizing the single particle Hamiltonian. 
We then calculated and stored in memory all the matrix elements H^' ^ = 



n = at + bU 



(30) 




These matrix elements can be constructed from: 



For la > 



n 



ai > 



(31) 



i=l 





1 



a-, ai nearest neighbors 
otherwise 



(32) 



(33) 
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These were then used to approximate the density matrix elements as: 

< a'l exp(-ei/)|a > ^ < a \{1 - eH)\a > (34) 

P was chosen large enough so that the results were converged to the exact 
result. For larger basis sets, the storage of matrix elements will not be 
feasible, but can be calculated "on-the-fiy" . During the course of the Monte 
Carlo calculation, we kept track of the N-electron states that appeared at 
each step. For the first 5000 steps, we periodically performed diagonalizations 
using only the states that appeared at least 0.1% of the time. Thus, a subset 
of the total number of N-electron states were diagonalized, using the current 
set of matrix elements H^' ^. After each diagonalization, all the Hamiltonian 
matrix elements were updated. Therefore, after a number of diagonalizations, 
the states we labeled as a were linear combinations of the original set of 
N-electron states. During subsequent Monte Carlo steps, we expected the 
adapted states to produce fewer kinks and to reduce the sign problem. 

This adaptive procedure is similar to stochastic diagonalization meth- 
ods,ii used to find the lowest eigenvalues of large Hamiltonians. The major 
difference is that here we include states that couple to excited states in the 
set of states used for diagonalization (which we must, since we are not focus- 
ing solely on obtaining ground state energies, but rather the partition func- 
tion). Stochastic diagonalization corresponds, roughly, to using the present 
approach, but limiting the number of kinks to 2 (one of the two states re- 
stricted to being the ground state) and performing a diagonalization every 
time a second kink is accepted. In addition, since the overall goal is to si- 
multaneously sample atomic positions and treat the electronic problem, the 
Hamiltonian matrix elements we are concerned with will change during the 
simulation, which is not the case in standard stochastic diagonalization tech- 
niques. In an application of the path integral approach to systems other than 
model systems, it is likely that many of the ideas of stochastic diagonalization 
can be "borrowed" for use in the adaptive scheme, while at the same time 
utilizing the power of the path integral method to treat finite temperatures 
and changing atomic positions. In any event, the adaptive procedure resulted 
in a new set of correlated N-electron states that were better representations 
of the eigenstates of the Hamiltonian. After diagonalization, the system was 
put in the ground state and the number of kinks was set to 0. After 5000 
steps, 15000 additional steps with no further diagonalization were taken in 
order to calculate the total energy. 
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The results are shown in Table 0. We first rewrote Eqn. ^ as 



QiP) = E(flE)pn(K}) (35) 

n=l \i=l ji I 

Here, n = 1 corresponds to the zero kink term and p is the density matrix 
corresponding to n and a specific set of states {a) We evaluated the average 
sign: 

,S) ^ <^ (36) 

For a poorly chosen set of states, (S*) will be nearly zero.i Thus, (S*) is a 
good measure of the severity of the sign problem. We also evaluated the 
average energy via 

/p \ '^-9pn/^|3) , . 

[-l^MonteCarlo) — ; ^ [O I ) 

\Pn) 

by taking the appropriate derivative of Eqn. An examination of Table ^ 
shows that the average sign of the density matrix is unity and demonstrates 
that the adaptive approach does indeed overcome the sign problem. An 
analysis of the number of types of kinks encountered indicated that only 
and 2 kinks were present after the 5000 step adaptive period. Since the 
density matrix is always positive for or 2 kinks, this lead to the value of 
1.0 for the average sign. Without the adaptive procedure, the average sign 
would have been very small and the usual sign problems would have been 
encountered. 

It is interesting to examine the convergence of the ground state energy 
to its final value during the adaptive phase of the simulation (the first 5000 
Monte Carlo passes), as well as the number of states involved in the diago- 
nalization. This information is given in Table |] and Table |^. It is evident 
that only a subset of the total number of states need to be included at any 
one diagonalization, which aids in the speed of the calculation. In addition, 
the energy is seen to converge relatively quickly to a final value during the 
adaptive stage. Neiectron = 3 and 5 doubly degenerate ground states, leading 
to 2 states being "visited" a substantial amount of time. Also note that, 
after the final diagonalization, the ground state energy of Neiectron = 4 is 
not the exact ground state energy, but the exact average energy is obtained 
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during the simulation with the addition of kinks. Finally, we note that the 
last diagonalization for N electron = 5, makes a very small change in the total 
energy. 

It may be noticed that, for each value of N electrons a large fraction of the 
states are involved in one of the diagonalization steps. This leads to the 
question of whether this is a required feature of the adaptive scheme. To test 
this, we performed an additional simulation for Neiectron = 5 and restricted 
the number of states that could be diagonalized at any time to 300. The 
results are shown in Table ^ from which it is evident that a diagonalization 
step involving a large fraction of the states is not needed for the adaptive 
scheme to work, although diagonalizing fewer states at each step does increase 
the length of the adaptive period. 

IV. Conclusions 

The adaptive approach has been shown to be successful in overcoming 
the sign problem inherent in standard path integral approaches. The use of 
the closed form expression for the kink expansion allows large values of P to 
be used. This method needs to be applied to additional systems, to assess 
its overall utility. In particular, there needs to be an assessment of how the 
method scales with the number of basis functions. This is not straightfor- 
ward to determine, since this will depend on what types of schemes are used 
to speed up the calculation. However, it is likely that, for fixed atomic posi- 
tions, that the method will scale in time similar to CI calculations, though it 
will not have to exhaustively evaluate all excitations (singles, doubles, etc.) 
and may therefore require less time than a traditional CI calculation. In 
addition, the size of the space of Slater determinants may grow too large for 
brute force diagonalization, while the adaptive approach can still be used 
(much in the same way that stochastic diagonalization can be used). Since 
the electronic calculation can be carried out simultaneously with atomic po- 
sition optimization, the time for a calculation with varying atomic positions 
should be comparable to the time required to perform a single point ah initio 
CI calculation. Additionally, as the adaptive wavefunctions become better 
representations of the exact wavefunction, the calculation should speed up, 
since only a small number of states will be coupled to the ground state. 
Future work will assess the applicability of this method to such types of 
calculations. 



15 



V. Acknowledgments 



Calculations were performed on computers purchased with a grant from 
the Louisiana Education Quality Support Fund. 



16 



References 

^C. Mak, R. Egger, and H. Weber-Gottschick, Phys Rev Lett 81, 4533 

(1998) . 

^C. Mak and R. Egger, J Chem Phys 110, 12 (1999). 

^R. Egger, L. Muhlbacher, and C. Mak, Phys Rev E 61, 5961 (2000). 

^S. Miura and S. Okazaki, J Chem Phys 112, 10116 (2000). 

^B. Mihtzcr, W. Magro, and D. Ceperley, Contributions to plasma physics 
39, 151 (1999). 

^W. Newman and A. Kuki, J Chem Phys 96, 1409 (1992). 

^R. Hall, J Chem Phys 94, 1312 (1991). 

^R. Hall and M. Prince, J Chem Phys 95, 5999 (1991). 

^D. Ceperley, Phys Rev Lett 69, 331 (1992). 
^°P. Roy, S. Jang, and G. Voth, J Chem Phys 111, 5303 (1999). 
^ip. Roy and G. Voth, J Chem Phys 110, 3647 (1999). 

^^B. Militzer, W. Magro, and D. Ceperley, Contrib Plasma Phys 89, 151 

(1999) . 

^^D. Ceperley, Path integral monte carlo methods for fermions, in Monte 
Carlo and molecular dynamics of condensed matter systems, edited by 
K. Binder and G. Ciccotti, 1996. 

^^R. P. Feynman and A. Hibbs, Quantum Mechanics and Path Integrals 
(McGraw-Hill, 1965). 

^^N. Rom, E. Fattal, A. K. Gupta, E. A. Carter, and D. Neuhauser, J. Chem. 
Phys. 109, 8241 (1998). 

^^N. Rom, D. M. Charutz, and D. Neuhauser, Chem. Phys. Lett. 270, 382 
(1997). 

i^Y. Asai, Phys. Rev. B 62, 10674 (2000). 

17 



'R. Baer, M. Head-Gordon, and D. Neuhauser, J. Chem. Phys. 109, 6219 
(1998). 

'J. Halting, O. Miiller, and P. Borrmann, Phys. Rev. B 62, 10207 (2000). 

'W. Hehre, L. Radom, P. Schleyer, and J. Pople, ab initio molecular orbital 
theory (Wiley 1986). 

K. Schotte and V. Schotte, Phys Rev B 4, 2228 (1971). 

R. Chiles, G. Jongeward, M. Bolton, and P. Wolynes, J Chem Phys 81, 
2039 (1986). 

I. Gradshteyn, I. Ryzhik, A. Jeffrey, and D. Zwillinger, Table of Integrals, 
Series, and Products (Academic Press, 2000). 

M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998). 

H. De Raedt and M. Frick, Physics Reports 231, 107 (1993). 



18 



Table Captions 



Tabic 1 Caption. Parameters for 3x3 Hubbard Model. Neiectron is the to- 
tal number of electron, Nup and Ndown the number of up- and down-spin 
electrons, respectively, (5 is the inverse temperature, < Eexact > is the exact 
thermally averaged energy, A is the difference between the ground and first 
excited states, P is the number of discretization points, < EMonteCario > is 
the average energy from the simulation, and < S* > is the average sign of the 
density matrix. Numbers in parenthesis represent 2 standard deviations. 
Table 2 Caption. Number of states included in the diagonalization for the 
adaptive steps performed during the first 5000 Monte Carlo passes. Diago- 
nalizations were performed after every 500 Monte Carlo passes up to 5000 
passes. In addition, diagonalizations were performed after 1 and 100 Monte 
Carlo steps. The total number of states are 324 (N=3), 1296 (N=4), and 
3024 (N=5). 

Table 3 Caption. Convergence of the ground state energy with adaptive 
diagonalizations. 

Table 4 Caption. Convergence of the ground state energy with adaptive diag- 
onalization for Nciectron = 5 and with the maximum number of states allowed 
to be diagonalized restricted to 300. After 5000 passes, diagonalizations were 
only performed when more than 2 kinks occurred for a significant amount of 
time. Following the adaptive period, the average energy was calculated and 
found to be -8.214186 (6). 
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-8.21418803 




42.6 


46.2 


36.0 
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524288 


< EMonteCarlo > 


-6.50920068(8) 


-7.649676(5) 


-8.21418801(4) 


< 5 > 


1.00(0) 


1.00(0) 


1.00(0) 



Table 1: 



Table |l] Parameters for 3x3 Hubbard Model. Neiectron is the total num- 
ber of electrons, N^p and N^own the number of up- and down-spin electrons, 
respectively, (3 is the inverse temperature, < E^xact > is the exact thermally 
averaged energy, AE is the difference between the ground and first excited 
states, P is the number of discretization points, < EMonteCario > is the aver- 
age energy from the simulation, and < 5* > is the average sign of the density 
matrix. Numbers in parenthesis represent 2 standard deviations. 
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Table 2: 



Table |. Number of states included in the diagonalization for the adaptive 
steps performed during the first 5000 Monte Carlo passes. Diagonahzations 
were performed after every 500 Monte Carlo passes up to 5000 passes. In 
addition, diagonahzations were performed after 1 and 100 Monte Carlo steps. 
The total number of states are 324 (N=3), 1296 (N=4), and 3024 (N=5). 
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Table 3: 



Table Convergence of the ground state energy with adaptive diagonal- 
izations. 
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Table ^. Convergence of the ground state energy with adaptive diagonal- 
ization for Neiedron = 5 and with the maximum number of states allowed to 
be diagonalized restricted to 300. After 5000 passes, diagonalizations were 
only performed when more than 2 kinks occurred for a significant amount of 
time. Following the adaptive period, the average energy was calculated and 
found to be -8.214186 (6). 
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Pass Number 


States Involved 


Energy 


1 


1 


-7.21199493 


100 


178 


-8.12195973 


500 


295 


-8.17901325 


1000 


257 


-8.18986053 


1500 


274 


-8.19522379 


2000 


275 


-8.19876874 


2500 


227 


-8.20054427 


3000 


283 


-8.20333449 


3500 


283 
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1000 


215 
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4500 


236 
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7000 
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6 
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3 


-8.21407828 
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-8.21412725 
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-8.21414318 
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3 
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4 
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Table 4: 
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